Interlayer magnetoresistance in multilayer Dirac electron systems: motion and 

merging of Dirac cones 
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We theoretically study the effect of the motion and the merging of Dirac cone on the interlayer 
magnetoresistance in multilayer graphene like systems. This merging, which could be induced by 
a uniaxial strain, gives rise in monolayer Dirac electron system to a topological transition from 
a semi-metallic phase to an insulating phase where Dirac points disappear. Based on a universal 
Hamiltonian proposed to describe the motion and the merging of Dirac points in two dimensional 
Dirac electron crystals, we calculate the interlayer conductivity of a stack of deformed graphene like 
layers using Kubo formula in the quantum limit where only the contribution of the n = Landau 
level is relevant. A crossover from a negative to a positive interlayer magnetoresistance is found 
to take place as the merging is approached. This sign change of the magnetoresistance could also 
result from a coupling between the Dirac valleys which is enhanced as the magnetic field amplitude 
increases. Our results may describe the behavior of the magnetotransport in the organic conductor 
a-(BEDT)2l3 at high pressure where the merging of Dirac cones could be observed. 

PACS numbers: 72.80.Vp, 72.15.Gd, 73.43.Qt, 71.10.Pm 



I. INTRODUCTION 

Since the discovery of graphene 1,2 , systems show- 
ing masslcss Dirac electron like dispersion relation con- 
tinue to attract considerable interest. The signature of 
such electrons has been recently brought out in the or- 
ganic conductor 3 a-(BEDT) 2 l3 where BEDT stand for 
bis(ethylenedithio)-tetrathiafulvalene. This compound 
consists of a stack of conducting BEDT layers separated 
by the insulating iodine planes. The weak coupling be- 
tween the conducting layers gives rise to the 2D character 
for the electronic properties of this material. Theoreti- 
cal studies 3,4 and band energy calculations 5,6 have given 
evidence for the presence of tow tilted Dirac cones which 
move under pressure. It has been argued that the in- 
terlayer magnetoresistance is a powerfull tool to probe 
the properties of the Dirac cones 7 . Tajima et al. 8 have 
observed a large negative interlayer magnetoresistance in 
a-(BEDT) 2 I.3 for a transverse magnetic field. This effect 
was ascribed to the carriers of the zero mode Landau 
level (n = 0). The authors have also reported a change 
to a positive magnetoresistance which was assigned to 
Zeeman splitting of the n — Landau level. 

A theoretical interpretation of these experimental re- 
sults has been proposed by Osada 9,10 by calculating, 
within a quantum approach, the interlayer magnetoresis- 
tance in a system of stacked Dirac electron layers. Osada 
showed that in the quantum limit, the negative magne- 
toresistance is due to the degeneracy of the zero mode 
Landau level which dominates the interlayer transport. 
At high field, spin splitting becomes relevant and gives 
rise to the crossover from negative to a positive mag- 
netoresistance due to the reduction of carriers density. 
Osada has also given an explanation of the angle depen- 
dence of the interlayer resistance observed by Tajima et 
al. 8 and which could not be understood within a semi- 



classical description. 

Based on transport measurements, Monteverde et 
al. have, recently, argued that the conduction in a- 
(BEDT)2l3 could not only be ascribed to Dirac electrons. 
Both massive and Dirac carriers contribute to the con- 
duction properties. Moreover, no merging of Dirac cones 
has been observed up to 3 GPa 11 . 

According to theoretical calculations 12 , the merging 
in a-(BEDT) 2 l3 is expected to occur around 0.5 GPa. 
More recently, Picchon et a/. 13 discussed, based on an 
analytical approach, the stability of Dirac points and the 
merging conditions in a-(BEDT) 2 I 3 . 

The motion and the merging of Dirac points has been 
observed in a tunable honeycomb optical lattice of ul- 
tracold Fermi gas 14 . Lim et al. 15 provided a theoretical 
description of the experiment of Dirac point manipula- 
tion in optical lattices. 

The topological phase transition from a zero gap band 
state, with two Dirac points, to a gaped phase was 
also observed in a microwave experiment simulating a 
strained graphene 16 . 

Montambaux et al. 17 have proposed a universal Hamil- 
tonian to describe the motion and the merging of Dirac 
cones in 2D systems. The proposed Hamiltonian of- 
fers the possibility the follow continuously the toplogical 
transition from the semi-metallic phase, with two Dirac 
points, to the insulating phase where the Dirac points 
merge and a gap opens in the energy spectrum. 

It is worth to stress that the merging in a monolayer 
graphene could not practically be observed since a very 
large strain is required 17,18 . It is then interesting to 
investigate the possible signature of this merging in 
other Dirac electron systems, in particular in multilayer 
graphene like system such as the organic conductor 
a-(BEDT) 2 I 3 . 



In this paper, we propose to study the merging of Dirac 
cones in a stack of undoped Dirac electron layers weakly 
coupled by a vertical tunneling. We look for the evidence 
of the merging in the behavior of the interlayer magne- 
toresistance. In the present work, we consider that each 
layer is described by the universal Hamiltonian proposed 
by Montambaux et al. and we introduce the interlayer 
coupling pcrturbativcly. For simplicity, we do not con- 
sider the tilt of Dirac cones which has been addressed in 
Refs.7 and 19. Moreover, we neglect the Zeeman effect 
which is, already, found to induce a sign change in the 
magnetoresistance 8 ' 9 . We also do not take into account 
the broadening of the Landau level which may result in a 
mixing of the Landau levels 9 . We derive, based on Kubo 
formula, the interlayer magnetoresistance in the quantum 
limit where only the zero mode (n = 0) is considered. In 
the next section, we focus on the behavior of the field and 
angle dependence of the interlayer magnetoresistance far 
from the merging of Dirac cones. In section III, we de- 
rive the magnetoresistance at the merging and discuss its 
experimental fingerprints. 



II. INTERLAYER MAGNETORESISTANCE: 
MOTION OF DIRAC CONES 

A. Independent Dirac valleys 



and move into the same direction to merge in a single 
point -Do denoted M in Fig. 2. 





We consider a stacking structure of layers weakly 
coupled along the transverse direction and we denote 
by t c the interlayer tunneling parameter (Fig. 1). Each 
layer is a grahene like system described by a triangular 
lattice with two atoms A and B by a unit cell. 



FIG. 2. (a) Deformation of the honeycomb lattice along the 
y direction. {0,1,0,2) is the lattice basis. The hopping pa- 
rameters to the first neighbors t and t' are different due the 
deformation, (b) Brillouin zone of undeformed graphene lat- 
tice. By applying a deformation, Dirac cones leave the K 
and K' points and move in the same direction and eventually 
merge in M point 18 . 




To describe the motion and the merging of Dirac cones 
in zero magnetic field, Montambaux et al. 17 have pro- 
posed the following Hamiltonian, so called, universal 
Hamiltonian 



H (P) = 




A + 2^ + ic vPv 



JLl. 

2)11 ' 



ICyPy 



(1) 



FIG. 1. Schematic representation of a multilayer system with 
magnetic field and current configuration, c denotes the inter- 
layer distance and 8 is the out of plane angle of the magnetic 
field from the conducting layer. 

In the absence of lattice deformation, the Dirac points 
D and D' = —D are in the points K and K' at the 
corners of the first Brillouin zone (BZ) 17 (Fig. 2). By 
applying a uniaxial strain, along the y direction for 
example, the Dirac points leave the corners of the BZ 



where p= (p x ,p y ) is the momentum measured relatively 
to the merging point Dq, c y is the electron velocity 
along the y direction, m* is an effective mass supposed 
to be positive and A is the parameter governing the 
topological transition. This two band Hamiltonian is 
written in the basis of the A and B site atom eigenstates 
(ip A ,ipB)- 

The universal Hamiltonian of Eq.l describes a de- 
formed graphene sheet in the presence of a uniaxial strain 
applied along the y axis 18 . The corresponding energy 



spectrum is: 




Plcl 



(2) 



Equation 2 shows an hybrid character, called by 
Montambaux et al. 17 a scmi-Dirac spectrum, with a 
Schrodinger like behavior along the x direction and a 
Dirac structure along the y axis. 

The case of A < corresponds to two distinct Dirac 
points along the x axis at ±pd where pr> = v— 2m*A 
whereas for A = 0, the Dirac points merge at Do. 

For A > 0, a gap of 2A opens in the energy spectrum 
and the system becomes insulator. 

To recover the full Dirac spectrum of undeformed 
graphene, the spectrum given by Eq.2 can be linearized 
along the x axis in the vicinity of Dirac points 17 . 

Let us now focus on the case where the deformed sys- 
tem is in the presence of a magnetic field 

B(B X = B cos cos <fi, By = B cos 8 sin <j>,B z = B sin 6) 



. The Hamiltonian given by Eq.l can be written, using 
the Peierls substitution p — > p + eA, as: 



#(7?) = 



2m* 



lCy7Ty 
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(3) 



where the effective momentum is given by if = (jr x = 
p x +ezB y —eyB z , ir y = p y —ezB x , 0) within the gauge A = 
(zB y —yB z , —zB Xl 0). n x and ir y satisfy the commutation 
relation [ir Xl Tr y ] = —iehB z . 

To derive the eigenfunctions and the energy spectrum 
of the Hamiltonian given by Eq.3, we consider for sim- 
plicity, as in Ref.17, the squared Hamiltonian H e ff. The 
eigenproblem reduces to: 



H, 



ipA 



= Ei 



ipA 

iJ'B 



(4) 



which may be written as: 
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(5) 
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where s = ± corresponds respectively to the A and the 
B sites. Equation 5 takes the form: 

[eg {P 2 y + e 2 z 2 B 2 x - 2ezB x P y ) + V(Y)] ^ B = E^l®, 

The potential V(Y) is given by 

,,he 2 B 2 



V(Y > = (0) <* + n 2 + .A 



rrr 



-Y (7) 



Here Y = syo — y, yo — Px e g" - and we introduce, as in 
Ref. 17 , the parameter 6 = Jr^-A. 

V(Y) has two minima at Y = ±\/\S\ separated by 
a distance of 2yfj[ which decreases as the transverse 
component B z of the magnetic field increases. 

V(Y) can be described by a two independent well for 
large \5\ corresponding to the case where Dirac cones are 
far from the merging point. In this case, an expansion of 
V(Y) around Fq yields to 



V(Y) - V(u) 



e 2 B 2 
2m* 



\5\u 2 ' 



s^^VWl 



here u = Y — Y and |u| <C \/\$\- The eigenvalues of Eq.6 
are then: 



E n 



± 



\/2hCy 



(8) 



hm* c v 



where I 

l\b\&Bl 

E n can be written in the form 1 ' E n = ±y / 2hc y c x eB z n 

as found by Himura et al. 19 in a-(BEDT)2l3 in the case 

of non tilted Dirac cones. 

In the limit of large negative 5, and using the notation 
of Ref. 9 , the eigenstates of Eq.6 corresponding to the 
layer position Zi take the following form : 

For the zero mode (n = 0) 



Fo,y ,zi(r) 







"0,1/0, Zi \i ) 



(9) 



For n > 0: 



±n,y ,z 



= 1 / ±h i Jf) \ 
y/2 V h n,yoMr) J 



where 



1 '7l,y ,Zi\' r ) /rr eX P 
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V exp 
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u„(y-y )5 z „ 



(11) 



Here L is the layer length along the x direction, yo is 
the center coordinate of the harmonic oscillator and 
u n (y — jjo) is the corresponding cigcnfunction. 

We now introduce the interlayer hopping Hamiltonian 
I 



AH 



-2i c cos^ as a perturbation 9 , c denotes the 



interlayer distance. 

From Kubo formula, and to the lowest order in i c , the 
interlayer conductivity a zz is given by: 



<r*«(w) 



ih 
V 



E E E 

spin n' ^'^z^ n,y a ,z, 



f(E n ,)-f(E n )\(n',y' ,zi\J z \n,yo,Zi}? 



E n > — E n 



fibj + ihr + E„ — E n i 



(12) 
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where J z is the interlayer current density J z = ¥ [z, AH] 
and t is the relaxation time. 

In the quantum limit, the mixing of Landau levels 
could be neglected and the relevant contribution to the 

I 



interlayer conduction is due to the zero mode Landau 
level n — 0. The matrix element of the J z can be ex- 
pressed in term of the effective interlayer hopping pa- 
rameter t r 9 : 



(o,y'o> z i\ J z\o,yo,Zi) = -r-tc 



,+ C 5~, ~ . By + S 



Z -J>~1 - 
' y Q ,yo- 



(13) 



where 



t c = t c cxp 



.eB x , 



y'o + yo 



cxp 



ec 
AHB~, 



-2-Bi + — Bt 



(14) 



The interlayer DC conductivity a zz = 5ft (o zz (uj = 0)) can 
then be written as 9 : 



2C 



t 2 c e 3 cr\B z 
irh 3 



-exp 



1 ec 2 



2hB z 
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2 D 2 

a B y 
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where the factor C given in Ref. 9 could be considered as 
a constant as far as the relaxation r is assumed to be 
held independent 9 . 
We denote by a the parameter measuring the amplitude 

of the strain: a - 



— =■ . The same parame- 

ter has been introduced by Himura et al. 19 as a measure 
of the the anisotropy strength. In our calculation, based 
on the universal Hamiltonian, this parameter describes 
the proximity of Dirac cones to the merging point. The 
smaller a, the closer the merging. 

In Eq.15, the proportionality of the prcfactor to B z is 
due to the Landau level degeneracy ^p- = |ff ■ 

Equation 15 is similar to the interlayer conductivity 
expression obtained by Osada 9 if one takes the case of 
undeformed graphene like system a = 1. 

It is worth to note that the expression of a zz given 
by Eq.15 is also reminiscent of that obtained by Himura 
et al. 19 in a-(BEDT)2l3 in the case of non tilted Dirac 



cones. 



If the hopping integrals are limited to the first neigh- 
boring atoms, the universal Hamiltonian parameters, in 
the case of fixed A and c v , take the form 17 : 



v —i - C"t» — 



v/3' 



«o 



t' 2 



andi 



3i'ao 
~2fr 



U6) 



where t' is the hopping parameter along the y direction 
which is different from the other first neighbor hopping 
parameter t regarding the effect of the uniaxial strain 
(Fig. 2). ao is the distance between the two atoms of the 
unit cell. The dependence of a on the hopping parameter 
t' is represented in Fig. 3. At the merging, t' = 2t, a 
vanishes. 

To derive the interlayer conductivity given by Eq.15, 
we have assumed a large negative A so that the two 
valleys of the double well potential V(Y) (Eq.7) could 
be considered as independent. The two Dirac cones are 
far from the merging point Dq. This assumption could 
be justified as far as a is not close to zero. One needs 
to define a criterion to fix the critical value of a below 
which the two potential valleys start to interact. This 
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FIG. 3. Dependence of the anisotropy parameter a = . ^- as 

a function of the hopping parameter ratio y . At the merging, 
a vanishes. 



FIG. 5. Interlayer conductivity as a function of the inplane 
component of the magnetic field expressed in the unit of -^ . 



point will be discussed later. 

According to Eq.15, a zz increases linearly with the field 
amplitude for a normal magnetic field as found in the 
undeformed case 9 . 

In the presence of an in-plane field component -By , 
the increase of a zz is reduced compared to the case of 
undeformed layer as shown in Fig. 4. This decrease is 
due to the Gaussian decay of the effective interlayer 
tunneling amplitude t c (Eq.14). The inplane component 
By generates a positive magnetoresistance effect since it 
induces an inplane Lorentz force. The latter reduces the 
interlayer tunneling giving rise to a positive magnetore- 
sistance. 
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FIG. 4. Interlayer conductivity as a function of the normal 
component of the magnetic field at a fixed inplane field com- 
ponent. 



sistivity p zz as a function on the magnetic field amplitude 
B. p zz is given by ' ' 



A 



-l^-A^Bl+a-Bl) 



(17) 



B + |-B z |exp 
Here Bq — 0.1 T is a fitting parameter 9 and A = 2 ct' i e: i c ' 
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FIG. 6. Interlayer resistivity p zz as a function of the field 
amplitude for different value of the anisotropy parameter a = 



' — . Calculations are done for 8 = 4 and 

c v 4 



? and the 



filed is expressed in the unit of 



Figure 6 shows a negative magnetoresistance for a 
close to unity. It turns out that the negative magnetore- 
sistance feature predicted in multilayer Dirac electron 
systems 9,19 is a robust effect which survives under 
uniaxial strain. 



According to Fig. 5, the decrease of magnetoresistance 
withBii is more pronounced as a decreases. 

In figure 6, we plot the dependence of the interlayer re- 



However, as a decreases, the Gaussian exponential 
decay in a zz (Eq.15) overcomes the linear increase with 
B z and a crossover from negative to positive interlayer 



magnctorcsistan.ee takes place at a critical field value 
B cr (a). The latter is found to decrease as a is reduced. 
This sign change of the magnctoresistance is due to a 
strain renormalization of the inplane field component Bh 
According to Eq.15, By is responsable of the Gaussian 
exponential decay giving rise to a positive magnctore- 
sistance. This effect is enhanced in the presence of 
the anisotropy which renormalizes the inplane field 
component: the effective field component, along the 
direction of Dirac point motion, B x is renormalized as 
^f which gets larger as the merging is approached (a 
decreases) . 

In figure 7 we present the dependence of the interlayer 
resistance on the field orientation from the conducting 
plane at fixed field amplitude and for different a val- 
ues. The maximum of the resistivity is, as in undeformed 
graphene layers 9 , for an inplane magnetic field (9 = 0). 
The results show that the peak for an inplane field is 
broadned as a decreases reflecting the enhancement of 
the inplane Lorentz force. This behavior is found to be 
unchanged by varying the azimuthal angle <f>. 






S 




— a=l 


0.33 




\ 


■■■■ a = 0.5 








- a = 0.3 


0.32 




\ 

71/4 ■ 

=5T 'v 


■- a = 0.2 


1 0.31 


6 = 
B 





~ 0.3 *- -. 



30 60 90 120 150 



♦ (deg) 



-20 20 40 60 SO 100 120 140 160 ISO 200 
8 (deg) 



FIG. 8. Dependence of the interlayer resistivity on the az- 
imuthal angle <j> at a fixed field amplitude for different values 
of a parameter. 



question. 

To summarize, the magnctoresistance shows, in the 
case of independent Dirac valleys, a <j) dependence for 
deformed honeycomb lattice and a possible sign change 
in the presence of an inplane field component. The 
latter gives rise to an inplane Lorentz force which gets 
enhanced as the strain amplitude increases (a decreases) . 
This enhancement results from a strain renormalization 
of the field component along the direction perpendicular 
to the deformation axis. 

At this point a natural question arises: below which 
value of a our assumption of independent Dirac valleys 
breaks down? In the following section, we try to bring 
some answers. 



FIG. 7. Out of plane angle dependence of the interlayer re- 
sistivity at a fixed field amplitude for different values of a 
parameter. 

The azimuthal angle dependence of the interlayer resis- 
tance for different values of a is depicted in Fig. 8 which 
shows that the motion of Dirac cones (a ^ l)gives rise to 
the $ dependence of p zz . The maximum of the resistivity 
corresponds to $ = where the inplane magnetic field 
component is aligned perpendicular to the strain direc- 
tion y. 

Morinari et al. ascribed the $ dependence of p zz in 
a-(BEDT)2l3 to the tilt of Dirac cones rather than to 
the anisotropy of Fermi velocity . We propose that a $ 
dependent interlayer magnetoresistance could be a signa- 
ture of Dirac cone motion. The question is how to dis- 
tinguish experimentally between the $ dependence due 
to the anisotropy of Fermi velocities (or the tilt of Dirac 
cones) and that induced by Dirac point motion? 

Fermi surface probes may be a key issue to unveil this 



B. Interacting Dirac valleys 

Regarding the double well structure of the potential 
ViY) (Eq.7), the Landau levels are doubly degenerate in 
the case where a is close to unity (large A). However, by 
approaching the merging (reducing a), the degeneracy is 
removed due to the tunneling between the Dirac valleys. 
The E n= o level splits into two levels separated by a gap. 
The system undergoes a crossover from a semi-metallic 
state, where E n= Q is degenerate, to a semi-conducting 
state where the degeneracy is lifted. This is expected to 
result in a change form a negative to a positive interlayer 
magnctoresistance since the carrier density of the lowest 
energy level is reduced. 

The valley degeneracy could also be removed, at 
a fixed value of A, by increasing the magnetic field 
amplitude: the distance between the potential minima 
(Eq.7) decreases with increasing the magnetic field 



/|A|, 



(2^0 ex 



,. A crossover from a negative to a positive 
magnctoresistance is expected at a critical field B z . cr (a) 
for a given value of the parameter a. 
As the merging is approached (a decreases), the two 
potential valleys are no more independent and a small 
field value could remove the valley degeneracy. The 
closer the merging, the smaller B z ^ cr . 

It is worth to stress that the crossover from the neg- 
ative to the positive interlayer magnetoresistance ob- 
served by Tajima et al. 8 in a-(BEDT)2l3 has been as- 
cribed by Osada 9 to Zeeman energy effect. Morinari and 
Tohyama 20 proposed that the sign change of the mag- 
netoresistance in a-(BEDT)2l3 is due to Landau level 
mixing effect which depends on the level broadening. 

In the following we argue that such crossover may 
be induced by a tunneling between Dirac valleys which 
takes place at high magnetic field or close to the merging. 

Montambaux et al. 17 have estimated the gap AE n be- 
tween the two energy levels obtained when the valley de- 
generacy of a Landau level E n is lifted. They found that 



-#i 



AE n 
Since a 2 = ^ = A /^PT, AE n takes the form 

c y Y m*c 2 ' ™ 

AE„ ~ c-^^e^c* 



(18) 



The expressions of m* and c y for fixed m* and A are 
I 



given by 



17 



III 



2h 2 

Sta 2 , 1 



3£a 



where a is the distance between neighboring atoms of 
the honeycomb lattice. 

AE n could then be written as: 



AE„ 



-4^2-5 



(19) 



where the dimensionless field is B z = 



B,t 



The critical field B ZiCr at which the gap AE n opens 
should scales, according to Eq.19, as: 



InB, 



6hu 



(20) 



It is interesting to note that the energy gap around 
the n = Landau level was also calculated by Esaki et 
al. who also found the same exponential behavior as 
Montambaux et al} 1 . 

To take into account the interaction between the Dirac 
valleys inducing the degeneracy lifting of the Landau 
level, we adopt a perturbative approach. We consider 
the first order correction to the Landau energy. The cor- 
responding wave function are those given in the previous 
section which are the zeroth order correction in terms of 
the perturbation 

The energy difference in the Kubo formula of the con- 
ductivity in the quantum limit (n = n' = 0) (Eq.12) is 
then replaced by AE n : 



(«) 






df i(o,&,z,qjs|o,j/o,zi)i !; 

dE hw + ihr + AE n 



(21) 
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The DC conductivity takes then the form: 



t\B, 



A 1 



t 2 AE 2 eXp 



1 ec l 



2hB z 



:Bi 



Q- 



^Bl 



where A 



2 J,% \ as given in the previous section. 

To have a rough estimation of the effect of the energy 
gap AE n on the conductivity, we make the following ap- 
proximation for the correction term 



-4V2^ 



1 



t 2 AE 2 
h 2 



where we replaced the energy gap AE n by its expression 
given by Eq.19. 

In figure 9 we plot the field dependence of the inter- 
layer conductivity a zz (Eq.22) in the case of a transverse 



magnetic field. The results show a change in the behav- 
ior of a zz at a critical field B z%cr which decreases with 
decreasing a. 

The dependence of B zcr on a is shown in Fig. 10 ac- 
cording to which lnB ZiCr ~ 6.3 In a which is in good 
agreement with the value estimated by Montambaux et 
al. 17 (Eq.20). This agreement support the approxima- 
tion of the first order energy correction we introduced in 
the Kubo formula (Eq.21) to account for the contribution 
of the Dirac valley tunneling to the magnetotransport. 

It turns out that for a given a, our assumption of in- 
dependent Dirac cones holds as far as B z < a 6 where 
B 7 



ih with B *° 



( 117. 



In a-(BEDT) 2 I 3 compound, B z0 - 660T since 
ao ~ 10 A which means that the assumption of inde- 
pendent Dirac cones is justified for a transverse field 
B, < 10T for a = 0.5. 
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FIG. 9. Dependence of the interlayer conductivity a zz on 
the transverse magnetic field B z for different values of the 
anisotropy parameter a. 



FIG. 11. Interlayer resistance p zz as a function of the field 
amplitude B for different value of the parameter a. The solid 
(broken) lines correspond to the results with interacting (in- 
dependent) Dirac valleys. 




FIG. 10. a dependence of the critical field B cr at which a 
change in the behavior of interlayer conductivity takes place 
(Fig. 9). The dashed line is a fit of the numerical data showing 
that \nB z cr ~ 6.3 In a. 



In Fig. 11 we plot the dependence of the magnetore- 
sistance p zz on the field amplitude B in the presence of 
an inplane field component. The solid lines correspond 
to the results obtained taking into account the coupling 
between the Dirac valleys whereas the broken lines are 
the results corresponding to the first section with inde- 
pendent Dirac cones. Fig. 11 shows that, in the presence 
of a coupling between Dirac valleys, the crossover from 
negative to positive magnetoresistance takes place at 
a smaller field compared to the case of independent 
Dirac cones. It turns out that the sign change of the 
magnetoresistance could be induced by Dirac valleys 
tunneling which is enhanced as the magnetic field 
increases (Eq.19) or as the merging is approached (a is 
reduced) . 

The outcomes of this section is that the crossover from 
negative to positive magnetoresistance could be induced 



by a coupling between Dirac valleys which gets more pro- 
nounced by approaching the merging or by increasing the 
magnetic field amplitude. 

In the following, we focus on the behavior of the inter- 
layer magnetoresistance at the merging and see how the 
field and the angle dependences of p zz are affected. 



III. INTERLAYER MAGNETORESISTANCE: 
MERGING OF DIRAC CONES 

The diagonalization of the universal Hamiltonian for 
A = and in the gauge A = (zB y — yB z ,—zB x ,0) re- 
duces to 

2 



HCy 

1 



IE 



Y 4 - 2sY 



i/>a,b = E^a,b (23) 



where we take, for simplicity, the squared Hamiltonian. 
The dimensionless operators H y and Y are given by: 

and 



n„ = i%, y 



r , Y 

7 ' 



yo - y, yo 



Px+ezBy 



■ B : 



ir y = py — ezB x . 7 is written as 17 
/ 2hc y m* 



7 



(24) 



Hy and Y satisfy the commutation relation 



Y, n, 



The eigenfunction of Eq.23 is of the form 
i/>(f) - c iciX c iC2y 0(Y) 



where a = *f* [y - z^j and c 2 = 2|*z. 0(F) is 
the eigenfunction of anharmonic quartic oscillator and 
y is the corresponding center coordinate. 4>(Y) is the 
eigenfunction of: 



flanh — 



he, 



n 



Y 4 - 2sY 



We will consider in the next the case of s = 1 since the Therefore, the solution 4>(Y) of the anharmonic part 

s = ±1 correspond to a symmetric problem as a function of Eq.23 can be written, for to the lowest Landau level 



of Y. 

In Refs.22 and 23, the authors studied the eigenprob- 
lem of the following anharmonic quartic oscillator 



and for Y > 0, as: 



<MF)=0o(f)~e-- 



(25) 



H ± = -—+g 2 x A ±2g\x\ 
The groundstate eigenfunction of the potential 

| x |3 

V- (x) = g 2 x A — 2g\x\ is of the form 4>q{x) ~ e" g_ 3~. 



The eigenstate of the zero mode level of Eq.23 takes 
then the form: 



Fo '^ (r1= U, io „(r1 







where 



(26) 



fo,y ,Zi(r) ~exp 



.eB z 



B„ 



l —\ y °- Zl B- ]X 



cxp 



i — z — y 



cxp 



(go - y) 

3 7 3 



r 



The interlayer hopping matrix is given, as in section II, by: 
I 



(Fo,y>z'\AH\F ,y a , Zi ) = -iciy^z'^yo^i) 



by> 7 -A- r b~., - . B v ~\- b 7 > 7 ._ r b~, - B v 



(27) 



r 



The effective interlayer hopping can then be written as: 
I 



i'.eB x (jjo+jk), , , 
tc(y >Zi,yo,Zi) ~i c exp( i— (Zi-Zi)) j dxexp 



■i i ,eB x ' 



exp 



1 ( x + a 



3 V 7 



exp 



1 ( x — a 



3 V 7 



with x = y — Vo ~L yo , a = Vo ~ yo = z *~ -g 2 - and 7 is given ductivity to the lowest order contribution of t 
by Eq.24 



o zz ~ti\B z \ 



cxp 



. eB x 

IC X 

n 



(29) 



Using the Kubo formula, we obtain the interlayer con- where c is the interlayer distance and f(x) is given by: 
I 



f(x) = exp 



1 / C \ fX By X 



3 V7/ \c 2B 



exp 



1 f C \ ( X By X 



3 V7/ \ c 2B 



(30) 



The normal field component B z appearing in the between the two atoms of the unit cell 17 . We then 
prefactor of a zz in Eq.29 is due to the degeneracy of the obtain M = f-B% where the dimensionless magnetic 
Landau level. 



/??'" 



In graphene like system, at the merging (A = 0), 

a where ao is the distance 



n l a and c lt 

3tai 



n 



field is B z = ^j-B z . We take for numerical calculations 
c = 1.75nm and a = IOA as in a(BEDT) 2 I 3 . 

The field dependence of the interlayer conductivity a zz 
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for a transverse magnetic field is plot in Fig. 12. 
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FIG. 12. Interlayer conductivity as a function of the nor- 
mal component of the magnetic field at the merging of Dirac 
points. 

Contrary to the case of separated Dirac cones, a zz de- 
creases, at the merging by increasing the normal field 
amplitude. This decrease appears as the continuity of 
the positive interlayer magnetoresistance obtained be- 
yond the crossover field B zcr for vanishing a in section 
II. We can then conclude that the crossover from nega- 
tive to positive magnetoresistance is due to Dirac cone 

motion. As a = ./— decreases, Dirac cones get closer 
y c v 

to the merging point and the interlayer conductivity is 

reduced. 

We do not claim that our calculations provide a contin- 
uous description of the Dirac point motion from the zero 
gap to the gaped phase. But one could notice that there 
is a sign change of the magnetoresistance, for a normal 
field {B = B z ), if the system moves from the indepen- 
dent Dirac valleys to the merging phase. In the former 
case, a zz is linear to B z whereas in the latter case a zz 
decreases with increasing B z . 

To obtain a continuous analytical description of the 
magnetotransport with the Dirac point motion, one need 
to derive the expression of the cigenfunctions of the an- 
harmonic quartic oscillator with the potential given by 
Eq.7. To the best of our knowledge, these eigenfunctions 
have not been analytically determined. Skala et al. 24 
found solutions for the Schrodinger equation correspond- 
ing to the Hamiltonian 



H 



-& + V ^ 



where V(x) — V1X+V2X 2 +V3X 3 +V4X 4 with the condition 
V4 > 0. The solution is of the form: 

ip(x) = exp (—gox — ffl^ /2 — g2% ) 

However, ip( x ) diverges for x — > ±00. 



In Fig. 13 we present the field dependence of the 
interlayer conductivity a zz in the presence of an inplanc 
field component. a zz shows a positive magnetoresistance 
as a function of the field amplitude. 



B(T) 



FIG. 13. Interlayer conductivity as a function of the magnetic 
field amplitude at the merging of Dirac points. 

Figures 14 and 15 show the field orientation depen- 
dence of the interlayer conductivity at the merging. 




30 
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FIG. 14. Dependence of the interlayer conductivity on the 
out of plane angle 9 at the merging of Dirac points. 

According to figure 15 the interlayer conductivity 
a zz shows a maximum along the strain direction y 
(3> = 7p. This behavior is found to be independent of 
the out-of-plane angle 9 as in the case where Dirac cones 
are far from the merging (Fig. 8). 

However, the dependence of a zz on the out of plane 
angle 9 is different from that found far from the merg- 
ing. The maximum of the conductivity is no more for 
a transverse magnetic field (9 = ?), as in Fig. 7, but is 
shifted towards 9 = as the inplane field component 
S11 is turned along the strain direction y. This behavior 
could be taken as an experimental probe for the merging 
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FIG. 15. Dependence of the interlayer conductivity on the 
azimuthal angle <j> for a fixed field amplitude at the merging 
of Dirac points. 



The latter induces an inplane Lorcntz force which reduces 
the interlayer tunneling. This effect is more pronounced 
as the amplitude of the deformation is increased. 

The motion of Dirac cones, resulting from the defor- 
mation, gives rise to a dependence of the interlayer mag- 
netoresistance on the azimuthal angle. However, the be- 
havior of the interlayer resistance with the out-of-planc 
angle 9 is unchanged compared to the undeformed case. 

We have argued that the sign change of the magnetore- 
sistance could also result from a coupling between Dirac 
valleys which removes the degeneracy of the Landau level 
and, hence, reduces the density of carriers. This coupling 
is enhanced as the merging is approached or at high mag- 
netic field. A criterion is proposed to define the range of 
validity of the independent Dirac valleys assumption. 

These features may be observed in a-(BEDT)2l3 under 
high pressure or in a stack of deformed graphenc like 
systems. 
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